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Abstract 

The principle of energy conservation leads to a generalized choice of transition probability in 
a piecewise adiabatic representation of quantum(-classical) dynamics. Significant improvement 
(almost an order of magnitude, depending on the parameters of the calculation) over previous 
schemes is achieved. Novel perspectives for theoretical calculations in coherent many-body systems 
are opened. 
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The importance of theoretical methods for the calculation of time-dependent quantum 
properties cannot be emphasized enough. The lack of general algorithms, so reliable as 
classical molecular dynamics simulations is to be contrasted with the manifold of open 
problems that scientists face both in condensed matter ^ and quantum information tech- 
nology jsj. Lately, we are also witnessing a renaissance of quantum approaches to biological 



phenomena a revival of interest generated by the combination of methodologies from 
open quantum systems |5] and quantum information theory ^ . Undoubtedly, the possibility 
of performing long-time quantum dynamical simulations would be an asset for all the above 
fields. 

When considering the calculation of time-dependent quantum properties, two main meth- 
ods are available: time- dependent density functional theory 0| and quantum-classical for- 
malisms Time-dependent density functional methods are usually limited to linear re- 
sponse while quantum-classical methods are restricted to perturbations around the adiabatic 
evolution, i.e., nonadiabatic corrections, of few relevant quantum degrees of freedom inter- 
acting with a classical bath. Nevertheless, quantum-classical methods promise to access the 

. Here, we are considering the 
formulation of quantum-classical theory by means of algebraic brackets which was proposed 
originally in and shown to arise from a linear approximation of the partially Wigner 
transformed quantum commutator [9]. It is remarkable that, when the environmental de- 
grees of freedom are harmonic and the coupling to the quantum subsystem is linear in the 
bath degrees of freedom, as in gauge theory [l^ , such theory becomes fully quantum. What 
is interesting from a computational point of view is that, within such a theory, a particular 
approximation (called momentum-jump in the adiabatic basis of the total system) leads to 
represent nonadiabatic dynamics in terms of piecewise (adiabatic) deterministic trajectories 



11 



121. 



interspersed by stochastic quantum transitions 

There is a long history of development and methods for treating non-adiabatic transitions 
with so called surface- hopping schemes Such schemes were originated in jisl. A more 
recent approach can be found in ^]. These methods are successful for the description of the 
dynamics but do not easily lead to an accurate formulation of the statistical mechanics of 
quantum-classical systems. Instead, the theory stemming from {q] allows one to address the 



consistent formulation of the quantum-classical statistical mechanics 



15| of general hybrid 



systems, i.e., the theory can describe, in the non relativistic limit, any quantum subsys- 
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tern coupled to a classical bath. It exactly conserves the energy of the total system and 
consistently describes the coupling between the quantum subsystem and the classical bath 
(or the quantum harmonic bath represented in Wigner phase space). The forms of the 
equations in the momentum-jump approximation also naturally provide a sampling tran- 
sition probability for nonadiabatic change of state. However, when nonadiabatic effects 
are included, the phase space trajectories representing the quantum(-classical) dynamics do 



no longer conserve the ener, 
time-step propagation 



Despite this, in its original formulation, called sequential 
12|, the algorithm is successful, although limited to somewhat 
short-times because of numerical instabilities arising from the sampling of the nonadiabatic 
transitions. The instability, in practice, restricts the range of applications of the method to 
charge transfers and rate processes 16j. 

More general quantum processes require the ability of sampling at longer time. In this 
letter we show how to achieve this by means of a suitable generalization (implementing the 
principle of energy conservation) of the transition probability in the sequential time-step 
propagation: this is the main theoretical idea we propose and it is introduced by Eq. (|5]) in 
the following. Before providing our solution, we sketch the theory and the original version of 
the sequential time-step propagation, which will be referred to in the following as primitive 
algorithm. The interested reader will find more details in 12|. The theory of quantum- 



classical dynamics is defined by the equation 
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where xw{X) is a quantum operator in a partial Wigner representation depending on the 
phase space point X = {R,P), comprised by coordinated and momenta respectively; Hw 
is the partially Wigner-transformed Hamiltonian operator of the total system, 13^ is the 
symplectic matrix, and d stands for the phase space gradient d/dX, with the arrow giving 
the direction of action of the operators. Without loss of generality, one can assume the form 
of the Hamiltonian to be HwiX) = ^ + hw{R) ■ In the adiabatic basis, defined by the 
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eigenvalue equation hw{R)\ci', R) = Ea{R)\a; R) , the quantum-classical evolution reads 

x'^"'{X,t) = J2{e''')^^,,,,x'/{X), (2) 

/3/3' 

where iCaa',i3i3' = iJ-'aa'^a/sSa'/s' + Jaa',pi3'- The diagonal operator iC^^^, = 
{ioJaa' + iLaa') Sa/s^a' 13' is defined in terms of the quantum adiabatic frequency Uaa'iR) = 
{Ea{R) — Ea'{R)) /h and of the classical-like Liouville operator iLaa' = (P/M) ■ d/dR + 
(1/2) {F^ + F^) ■ {d/dP), where F{y are the Hellmann-Feynman forces f2d]. The operator 
Jaa',i3ii' IS purcly off-diagoual and its action realizes the quantum nonadiabatic transitions. 
It is worth remarking that Eqs. ([T]) and ([2]) exactly conserve the total Hamiltonian of the 
system Hw{X). 

In the momentum-jump approximation, the form of the off-diagonal operator Jaa'^pp' is 
simplified 12]. Here, we denote such an approximation by J^^\pi- The action of J^f^/^/^/ 



1 momenta. The technical details can be 



changes the quantum state and rescales the bat 
found, among many other possible references, in 

one can also define a momentum-jump Liouville operator, iC^^^g, = iC^^' bb' + J^*-™^^ 



12| . Using the momentum-jump operator, 
(mp) — ^ro _|_ T(mp) 

aa',BB' ~ ''^aa',BB' ' "^aa^BB'^ 



approximating the exact operator iCaa',BB' ^'i- (ED- 

Deterministic dynamics is too-complicated to be solved, so one has to resort to stochastic 
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schemes. A very elegant approach is provided by the sequential short time propagation 
(the primitive algorithm). This is summarized as follows. For a small time step r the 
quantum-classical propagator is approximated as 



One can prove that the concatenation of short time steps, according to Eq. ([3]), reproduces 
exactly the Dyson integral expansion of the operator exp ^^/- The algorithm 

unfolds by considering the action of J(™p) dictated by a stochastic process with a certain 
transition probability. The form of J(™p) naturally suggests the following primitive choice 
of the transition probability (for example considering the a — > /3 quantum transition): 

I + ■ d^p{R)\T 

where daB = («; R\d/dR\(3] R) is the coupling vector. Normalization fixes the probability of 
not making any transition in the time interval r as Qi°j(X, At) = 1 - r^^l The stochastic 
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propagation amounts to deterministic trajectory segments, propagating on single or mean 
energy surfaces, interspersed by transitions between energy surfaces. The transitions break 
the conservation of energy along the single trajectory: the conservation is only satisfied in 
an averaged sense in the ensemble. As one can see from Eq. (jlj), arbitrary shifted momenta 
P' can arise from a sampled transition. As experience has shown, this leads in general to 
very big denominators in the left hand side of Eq. This denominators get multiplied 
with each other along the trajectory to give its overall weight in the stochastic ensemble. 
The concatenation of big weights arising from nonadiabatic transitions produces a numerical 
instability which has so far limited the application of the primitive algorithm to somewhat 
short times. 

The principle of energy conservation, which is exactly satisfied by Eq. ([1]), guides us in 
defining a generalized transition probability as 

1 + T\{a\(5)\w {c£,Saa',l3l3') 

where we have defined {P/M)-dap{R) = Again, normalization provides At) = 

1 — 'Paid- Upon introducing the variation of energy in any quantum transition Eaa'^pp' = 
^ + \ iEa{R) + Ea'{R)) ~ ^ ~ \ iE(s{R) + Epi{R)), the generalized weight introduced in 
Eq. ([5]) is defined as 

I 1 if \£aa',l3l3'\ < Cs] 

w [ce) = ( (6) 
I otherwise; 

with C£- tunable constants controlling the numerical error on the energy conservation. 

The generalized transition probability in Eq. (|5]) and the energy-conserving weight in 
Eq. ([6]) are our fundamental findings, improving the primitive algorithm. Because of our 
choice of w (c^:), the sampled transitions can only allow shifted momenta P' which conserve 
(within the required numerical error specified by ) the energy of the system. This in turn 
avoids arbitrarily big denominators in Eq. and dramatically improves the stability of the 
algorithm. 

In order to illustrate the efficiency of the energy-conserving sampling, we performed a 
series of calculation on the dynamics of the spin-boson model [2]]]. The partially Wigner 
transformed Hamiltonian of this model (in scaled coordinates) is Hw = —^^x+J2f=iiE'j / s+ 
co'jR'j/2 — CjRjCTzi where and are the Pauli matrices, Rj,Pj are coordinates and mo- 
menta of harmonic degrees of freedom (in the following we have used A^ = 200, Cj are 
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FIG. 1: 

Comparison of the primitive (A) and energy-conserving (•) sampling for f3 = 0.3, = 1/3, and 
^ = 0.007. The calculation with the primitive algorithm are propagated until t = 20 and then 
stopped since already around t = 15 the statistical error becomes very big, as the error bars show. 
The calculation with the energy-conserving sampling (•), with C£ = 0.01, can be extended further 

than t = 30. 

coupling constants defined in term of the Kondo parameter £). Details of the model and 
definition of coordinates and parameters can be found in 12|. Figures [1] and |2] illustrate 
the numerical comparison between the primitive and our energy-conserving sampling for the 
relaxation dynamics of for various couplings, temperatures, and tunnel splitting Q. The 
results obtained with the primitive algorithm are displayed by white triangles while those 
obtained with our energy-conserving sampling by black filled circles. Figure [1] displays the 
results of the numerical calculation for /3 = 0.3, Q = 1/3, and ^ = 0.007. Basically, after 
t = 15 (in dimensionless units) the error bars on the primitive algorithm results grow ex- 
ponentially fast and the calculation is stopped at t = 20. Instead, the calculation with our 
generalized sampling scheme can be extended further than t = 30: for this set of parameters 
we obtain an improvement over the time interval we can sample of at least two. Figure [2] 
displays the results of the numerical calculation for /3 = 12.5, Q = 0.4, and ^ = 0.09. This 
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time, the statistical errors of the primitive algorithm start growing fast around t = 10, while 
our scheme can reach further than t = 100, providing an improvement of an order of mag- 
nitude. Summarizing, the simulation shows that the use of our energy- conserving sampling 
dramatically improves the stability of the elegant sequential time step algorithm at long 
time. 

In conclusion, it is worth mentioning that the approach, embodied by Eq. to modify 
the transition probability in order to respect a conservation law and improving the stability 
is very general: it is by no means restricted to quantum(-classical) dynamics in the partial 
Wigner representation. On the contrary, there are reasonable expectations that the gener- 
alized scheme that we have presented here can be applied, after suitable changes, to other 
stochastic approaches for calculating time- dependent properties, both in the classical and 
quantum cases. 
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FIG. 2: 



Comparison of primitive (A) and energy-conserving (•) sampling for /3 = 12.5, = 0.4, and 
^ = 0.09. The calculation with the primitive algorithm are propagated until t = 20 and then 
stopped since already around t = 12 the statistical error becomes very big, as the error bars show. 

The calculation with the energy-conserving sampling (• and a continuous line to help the eye), 
with C£ = 0.01, can be extended further than t = 100 (the error bars of the order of magnitude of 

the • symbols). 
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